Macroscopic resonant tunneling of magnetic flux 
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We have developed a quantitative theory of resonant tunneling of magnetic flux between discrete 
macroscopically distinct quantum states in SQUID systems. The theory is based on the standard 
density-matrix approach. Its new elements include the discussion of the two different relaxation 
mechanisms that exist for the double-well potential, and description of the "photon-assisted" tun- 
neling driven by external rf radiation. It is shown that in the case of coherent flux dynamics, rf 
radiation should lead to splitting of the peaks of resonant flux tunneling, indicating that the resonant 
tunneling is a convenient tool for studying macroscopic quantum coherence of flux. 



I. INTRODUCTION 



Bose condensates of Cooper pairs in superconductors 
have a remakable ability to populate a single quantum 
state with a macroscopically large number of particles. 
This property of the Cooper pair condensates leads to 
macroscopic quantum effects in the dynamics of Joseph- 
son junctions, i.e. makes it possible for a Josephson junc- 
tion to behave as a pure quantum state of a simple quan- 
tum mechanical system while still containing macroscop- 
ically large number of Cooper pairs. Due to the strongly 
nonlinear character of the Josephson dynamics even in 
the classical regime, macroscopic quantum effects can be 
quite non-trivial. They have, been known for the past 
twenty years - see reviewsantl, and are continuing to at- 
tract considerable interesta. This interest is to a large 
extent stimulated by the challenge to our understanding 
of the foundations of the quantum theory presented by di- 
rect manifea|ta|tions of quantum mechanics at the macro- 
scopic leveflij. Another recent motivation for studying 
quantum effect in Josephson dynamics is prftvided by 
possible applications to quantum computationB'EI. 

The most advanced macroscopic quantum effect ob- 
served experimentally up to now is resonant tunneling 
between quantized energy-Iejvels in the adjacent wells of 
the Josephson potentialEJO. The aim of this work is 
to develop a theory of this phenomenon in SQUID sys- 
tems, where the potential contains two wells each with 
a different value of the average magnetic flux. We con- 
sider the regime of weak energy dissipation important 
for studying coherent effects in resonant tunneling. This 
regime has not been discussed appropriately in the ex- 
isting tresLtments of resonant tunneling in double-well 
potentialsE£rEj or multi-well potentials papcisponding to 
the current-biased Josephson junctionsEjO. The most 
essential new feature of our approach is an account of 
the two types of relaxation mechanisms, intrawell and 
interwell, that exist in the system. The two relaxation 
mechanisms are very different in their dependence on the 
parameters of the SQUID potential, and lead to different 
shapes of the resonant tunneling peaks. Differences in 
relaxation mechanisms also make macroscopic resonant 
tunneling of flux different from the otherwise very similar 



"mesoscopic" resonant tunnslmg. between charge states 
of small JoseuhscLo junctionsO'tS and electron states in 
quantum dotal3~Ej. 

Another new element of this work is the discussion 
of the "photon-assisted" macroscopic resonant flux tun- 
neling under rf irradiation. We show that in contrast 
to tunneling under stationary-bias conditions, the peaks 
of the photon-assisted tunneling depend qualitatively on 
the strength of decoherencc in the flux dynamics. In the 
case of coherent flux dynamics, the resonant peaks of the 
photon-assisted tunneling are split in two. The splitting 
reflects the coherent hybridization of the macroscopic 
flux states in the two wells of the SQUID potential and 
is suppressed with increasing rate of decoherence. Very 
recently, such a splitting of the resonani-Jlux-tunncling 
peaks has been observed experimentallycj, demonstrat- 
ing the quantuijn_|Coherence of the macroscopically dis- 
tinct flux statestSo. The paper is organized as follows. 
In Sec. 2 we derive the evolution equations for the den- 
sity matrix describing the resonant flux tunneling under 
stationary-bias conditions, and introduce the two relax- 
ation mechanisms for tunneling dynamics. Using these 
equations, we calculate the rate of flux tunneling in Sec. 
3. In Sec. 4, we extend the results of Sees. 2 and 3 to the 
case of photon-assisted tunneling. 

II. EQUATIONS FOR THE DENSITY MATRIX 



To derive equations for the density matrix in the 
regime of resonant tunneling of magnetic flux $ we 
consider the standard model of the phase dynamics in 
SQUIDS. The combination of the magnetic energy of the 
SQUID loop biased with an external flux and the Joseph- 
son coupling energy of the SQUID junctions (for details, 
see e.g.j2il) produces the double- well potential t/(<&) for 
$ evolution (shown schematically in Fig. 1 below). The 
main part of the Hamiltonian governing the flux dynam- 
ics consists of the potential U ($) and the charging energy 
of the junction capacitance C: 
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The charge Q on the junction capacitance and the flux 
$ satisfy standard commutation relations [$, Q] = ih. 

The two wells of the potential J7($) have discrete en- 
ergy states Ejn with characteristic energy separation on 
the order of Wj, where j = 1,2 is the well index, ojj are 
the oscillation frequencies around the potential minima, 
and n — 0,1,... numbers the states within each well. 
The two frequencies ujj have the same order of magni- 
tude LOi ^ L02 = ujp- External magnetic flux controls the 
energy difference between the states in opposite wells. 
Away from the resonance conditions, when all the ener- 
gies Ejn are separated by large energy gaps of order ojp, 
the states \ jn) are localized within the jth well, and the 
amplitude of the wavefunctions '!/'j„(4>) in the opposite 
well is very small. However, when the energies of the two 
states |1) = |lni) and |2) = |2n2) are close, \e\ <C Wp, 
where e = eim — £2712 7 these states become strongly cou- 
pled, and the wavefunctions spread over the both wells. 
As shown in the Appendix, strong coupling of the states 
|1, 2) at resonance can be described by the tunneling am- 
plitude A, and the Hamiltonian (||) reduces to the regular 
two-state form in the basis of these states: 

^o = i[e(|l>(l|-|2)(2|)-A(|l)(2|-|-|2)(l|)], (2) 

where D is the quantum mechanical transparency of the 
barrier separating the wells. 

Without perturbations, the two-state dynamics de- 
scribed by (||) is decoupled from the other states of the 
Hamiltonian (|^). The most important perturbation cre- 
ating such coupling is the energy dissipation that induces 
transitions between the states |1, 2) and other states \ jn). 
In the relevant temperature range below superconduct- 
ing energy gap of the junction electrodes, the quasiparti- 
cle tunneling is suppressed and the main source of energy 
dissipation is the electromagnetic environment of the sys- 
tem. Under the assumption that the electromagnetic 
modes of the environment are in equilibrium at tempera- 
ture T, and are well described by linear electrodynamics, 
the interaction between the flux $ and the heat bath of 
these modes can be written as 
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Here // is the fluctuating current created in the SQUID 
loop by the environment, with the correlation function 
given by the fluctuation-dissipation theorem: 



{If{t)If{t + T)) = 
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where the brackets (...) denote averaging over the equilib- 
rium density matrix of the environment, and G(w) is the 
dissipative part of the environment conductance. Equa- 
tion is sufficient to characterize completely the ef- 
fects of the weak energy dissipation considered in this 
work. For arbitrary dissipation strength, one can use 



the Caldeira-Leggett modeled to express explicitly the 
environment Hamiltonian and the current operator If in 
terms of a set of harmonic oscillators. 

The interaction (||) induces both "vertical" transitions 
within each well and direct intcrwcU transitions. In terms 
of the two-state dynamics with the Hamiltonian (||), the 
latter correspond to the modulation of the tunneling am- 
plitude A by the environment. The matrix elements of 
this type of interwell transitions are, however, smaller 
by a factor of A/up than those of the intrawell transi- 
tions. The small matrix elements can be neglected un- 
der the conditions of resonance, when the flux tunneling 
between the two wells is dominated by the stronger res- 
onant processes. In this approximation, we can omit the 
terms in the flux operator in the interaction (^) that are 
non-diagonal in the well index j: 



(5) 
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Here are the matrix elements of $ in the \ jn) ba- 

sis. The perturbation with the flux operator (j^) has 
two effects on the dynamics of the states |1, 2). The first 
is fluctuations of the energy difference e, which induce 
transitions between these states and lead to the loss of 
mutual coherence between them. The part of the $ oper- 
ator (^) responsible for these fluctuations can be written 
as (5$( 1)(1| - |2)(2|), where is half of the difference 
between the average flux values in the states |1) and |2). 
The remaining terms in (||) induce intrawell transitions 
from the states |1) and |2) to the other states in the cor- 
responding wells. 

For weak dissipation both effects can be described 
quantitati^ly by the standard density matrix technique 
- see, e.g.,E3. The description starts from the equation for 
the evolution of the density matrix p, obtained treating 
the coupling V (^) in second-order perturbation theory: 

m = -t[Ho, p] - 1* dT{[Vit), [V{t), p(t)]]) . (6) 

If the environment has a large cut-off frequency loc S> 
£, A, the density matrix evolves slowly on the time scale 
of variations of V{t), and we can make the Markov ap- 
proximation p(r) ~ p{t) in the last term of eq. The 
condition of weak dissipation also allows us to keep only 
the dissipative terms in eq. (^) that do not oscillate in 
time with frequencies of the main Hamiltonian Hq, since 
only these terms lead to effects that accumulate with 
time. Using these approximations to evaluate the dissi- 
pative part of the eq. (|^), we obtain the final equation for 
the density matrix in the basis of resonant states |1,2) 
relevant for the transfer of flux between the wells: 



p = -i[Ha,p] + T[p] +7[p] . 



(7) 



The term F in this equation describes the effect of the 
intrawell transitions from the states |1) and |2): 
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At temperatures smaller than energy separation in the 
wells, the total decay rate Tj in eq. (^) of the state \j) 
into the states with lower energy in the same well is: 

2|(j)0')^ |2 

= ^ — {Snj — £n)G{enj — £n) ■ 

n<nj 

The second dissipation term "f[p\ in eq. describes 
transitions and decoherence within the |1), |2) subspace, 
and has a simple form in the basis of energy eigenstates 
of the two-state Hamiltonian (||): 

7[P] = -C/t f ^^+"-^>'Au. (9) 

Here r is the density matrix in the eigenstate basis: R = 
U pW , and U is the rotation matrix from this basis to 
the flux basis |2): 

U=[il- el^Yl^o, + (1 + e/17)V2^^]/V2 , 

where the ct's denote Pauli matrices, and = (e^ + 
^2^1/2^ The transition rates ^d,u and the decoherence 
rate 7 are: 
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7u = 7de 



-O/T 



with the dimensionless parameter g = 2G((5$)^/?i char- 
acterizing the strength of the interwell relaxation, and 
we assumed that G[ui) is constant in the small-frequency 
range uj ^ SI. 

Equation (|^), with the dissipation terms (||) and (|), 
is used below to describe resonant tunneling of flux be- 
tween the two wells in various regimes. Before doing 
this, we discuss the relative magnitude of the two dis- 
sipation terms in this equation. Since the width of the 
wells is of the same order of magnitude as the barrier be- 
tween them, the magnitude of the flux matrix elements 
of the intrawell relaxation r[p], and of the interwell re- 
laxation 7[p] (determined, respectively, by the "width" 
of the wavefunctions inside the wells and the distance 
between the wells) should be close. The main difference 
between the two relaxation mechanisms is that the in- 
trawell transitions dissipate energy lOj , whereas the inter- 
well relaxation 7[p] involves only much smaller energies 
on the order of e. A, T . This means that the intrawell 
relaxation typically dominates the flux-tunneling dynam- 
ics. In particular, even under the assumed condition of 
weak relaxation (which for r[p] means that the rates V j 
are small compared to the oscillation frequencies uij) the 
rates Fj can still be much larger than the frequencies 
e, A of the two-state dynamics. The interwell relaxation 
7[p] will generally only be stronger than the intrawell 
relaxation if the environment has relatively low cut-off 
frequency ^ t^p. 



III. STATIONARY BIAS 

In this Section, we calculate the rate of the resonant 
tunneling between the wells in the situation when the 
external flux through the SQUID loop does not contain 
an ac component, and the energy dissipation drives the 
initial flux state in the left well towards equilibrium. At 
low temperatures T <C the flux stays in the ground 
state of the left well and tunnels into the right well out 
of this state (Fig. 1). In this case, the relaxation rate Fi 
obviously vanishes. We begin by considering the situa- 
tion with only the intrawell relaxation present. Equation 
(0) can then be written in matrix elements as: 

pii = AImpi2 , P22 = -AImpi2 - 2rp22 , 

P12 = -{%e + F)pi2 + (iA/2)(p22 - Pii) , (10) 

with r2/2 = r in this Section. 



U(^>) 




FIG. 1. Schematic diagram of the macroscopic resonant 
tunneling of flux <1> in the double-well potential {/('I?) at low 
temperatures under stationary-bias conditions. The flux tun- 
nels out of ground state |1) in the left well coupled with the 
amplitude —A/2 to the resonant state |2) in the right well, 
where it decays with the rate r2 into the lower states in this 
well. 

After transformation to the real and imaginary 
parts of the off-diagonal matrix element pi2 and the 
sum/difference pn ± P22 of the diagonal elements of the 
density matrix, eqs. ( p^ ) can be solved directly with the 
initial condition that the flux is in the left well at time 
i = 0, pii(O) = 1: 



-rt 



(r2+t^2^((l + — )coshAt-H 



2- sinh At) + (A^ - r2)((l - ^^) costjt - 2 - sinwt) 
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(11) 
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In eq. (|1 1|), the eigenfrequencies uj and A of the system 
of equations (noh are: 



cj , A = 



-,1/2 



(12) 



Equations (|ll| ) contain aU the information about dynam- 
ics of the flux tunnehng. When the relaxation rate V2 is 
much smaller than the oscillation frequency Q,, the tun- 
neling process consists of the weakly damped coherent 
oscillations of the flux between the wells followed by re- 
laxation in the right well. With increasing relaxation 
rate the oscillation part of this process becomes increas- 
ingly more damped and turns into incoherent jumps of 
the flux from the left into the right well represented by 
the non-oscillatory exponential decay of pn. 

The matrix element pii{t) has the meaning of the prob- 
ability for the flux to remain in the left well at time 
t. In the case of monotonous decay of pn, its deriva- 
tive f{t) = — [pii{t)]' gives the probability density of the 
time for flux tunneling between the wells. This proba- 
bility density can be used to find all statistical charac- 
teristics of the flux tunneling process. For instance, the 
average time r it takes for the flux to tunnel can be calcu- 
lated as r = dt tf{t) = Jp°° dtpii{t). In the coherent 
regime, when pii(i) oscillates in time, f(t) can be neg- 
ative and cannot be interpreted as probability density, 
and the question of how to define the tunneling time r 
becomes non-trivial. To define t, one needs to establish 
the event that terminates the tunneling process. The def- 
inition adopted below assumes that the tunneling process 
ends when the fiux makes the transition from the state |2) 
into one of the lower energy states in the right well. This 
definition is motivated by the fact that such a transition 
eliminates the possibility for the flux to return into the 
left well. With such definition, the time the flux spends 
in the state |2) is included in the tunneling time r which 
then should be calculated as 



dt{pn{t) + p22{t)) . 



(13) 



From equations (|11|) and (y_3|) we obtain the tunneling 
rate: 



2A2 



Tl + 4^2 • 



(14) 



Equation (|lj) describes the Lorentzian peak of the res- 
onant flux tunneling. It shows that the resonant flux 
tunneling under stationary-bias conditions does not al- 
low one to distinguish qualitatively between the regimes 
of coherent and incoherent flux tunneling since the shape 
of the resonance peak ( p^ remains the same regardless 
of the magnitude of the relaxation/decoherence rate F. 

The average tunneling rate can be calculated with- 
out explicit solution of the time-dependent equations for 
the density matrix. Instead of attempting to describe an 



individual tunneling event with the time-dependent so- 
lution, we can consider a large number of these events, 
assuming that after each transition from the left to the 
right well the system is immediately returned back to its 
initial state and the process is repeated. This immediate 
return means that the system is effectively decaying from 
the resonant state in the right well directly into the ini- 
tial state in the left well, and can be modeled by adding 
the term F2/922 into the equation for pu. With such a 
modiflcation, eq. (^ has a non-trivial stationary solution 
p^^\ and the tunneling rate r^^ defined by eq. (|l3|) can 
be found from this solution as 



-1 (0) 

r = T2py 



(15) 



Since this method does not require solution of the time- 
dependent equations for the density matrix, it consider- 
ably simplifies the calculation of the average tunneling 
rate. 

To illustrate this procedure, we consider first the same 
tunneling under the stationary bias conditions described 
by eqs. ( [lO| ) with the term F2P22 included into the equa- 
tion for pii. Solving the stationary equation for the off- 
diagonal element pi2 of the density matrix and plugging 
the solution into the equations for the diagonal elements, 
we get the simple rate equations: 



Pii = r'(P22 - /Oil) + F2/922 , P22 = - 



Pii 



(16) 



where the transfer rate F' = A2F2/(F^ 



. 2 -r 4£2) between 
the two wells can be viewed as the "Golden-rule" rate 
of transition with the matrix element A/2 into the state 
1 2) broadened by the relaxation F2. From the stationary 

solution of eq. ( pi) w e find that p^22 = r7(r2 + 2F'), 
and see that eq. (1^) indeed reproduces the tunneling 
rate (|l|). 

The tails of the resonant peak (|lj) at e ^ A, F2 allow 
for another simple interpretation. At large e, the wave- 
function of the eigenstate of the two-state Hamiltonian 
localized in the left well has the probability amplitude 
A/2e in the right well. The tunneling rate F in this 
regime can be found then as the probability to be in the 
right well times the relaxation rate F2: 



F = F2AV4e^ 



(17) 



This simple reasoning indeed reproduces the tails of the 
peak (|l^) and allows us to obtain an estimate of the 
tunneling rate between the resonances. As shown in the 
Appendix, the wavefunction amplitude in the right well 
between the resonances is 7rA/[2w2 sin(7re/w2)]. From 
this we can write: 



= F2 



ttA 



2w2 sin(7re/aj2) 



(18) 



It should be noted that eq. ( [l^ ) is only an estimate, since 
F2 and A depend on e for e ~ u;2, and can be different 
from their values at resonance. However, ai e ^ lu2, 
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they are constant, and eqs. ( p^ ) and ( [l8| ) coincide for 
e ^ A,r2. In this range of e the tunnchng rate r"^ 
changes as e^^. 

If the interwell relaxation is non-neghgible, we need to 
keep both relaxation terms in eq. (j^). The assumption 
that the relaxation rates are small in comparison with 
A, allows us to use eq. (^) for the interwell relaxation 
and makes it convenient to consider the flux dynamics in 
the eigenstates basis. The stationary values of the off- 
diagonal elements of the density matrix r in this basis 
are vanishing for weak relaxation. To find the diago- 
nal elements of r, we transform eq. (^) for the intrawell 
relaxation (with Fi = and added term T2P22 in the 
evolution of pn) into this basis. Neglecting rapidly oscil- 
lating terms, we see that the diagonal part of the weak 
intrawell relaxation is: 



rn=r.[^ + i(l + l,)(r.,-rn)] 



r22 



-rii 



Combining this expression with eq. (j^) for the interwell 
relaxation, we find the stationary values of the diagonal 
elements of r. After transformation back to the flux ba- 
sis we finally obtain the stationary element p^22 
density matrix p in the flux basis and the flux tunneling 
rate (|l|): 



17coth(f7/2T) 



P 



p = T2/2g. 



(19) 



The parameter p in eq. ([l9| ) can be interpreted as 
the energy at which the characteristic interwell relaxation 
rate (which increases with the energy difference between 
the two resonant states) becomes equal to the rate r2 of 
the relaxation in the right well. The tunneling rate ( |l^ ) is 
plotted in Fig. 2 for zero temperature and several values 
of p. The interwell relaxation becomes stronger with de- 
creasing p, making the the resonant-tunneling peak pro- 
gressively more asymmetric. For negative bias e, when 
the state |2) in the right well is higher in energy than the 
state |1) in the left well, the interwell relaxation is sup- 
pressed for low temperatures, T <C A, and the tail of the 
resonant peak (^9|) coincides with that of the Lorentzian 
peak (|^. On the other hand, for positive £, transitions 
from the state |1) into |2) are allowed, and if < A, 
the interwell relaxation dominates for all positive e. The 
tunneling rate decreases in this case only as 1/e with 
increasing e. 

At e 3> A, and vanishing temperature T, the tunneling 
rate ( [l9| ) determined by the interplay between the inter- 
well and intrawell relaxation can be understood in terms 
of the competition between the two tunneling paths (see 
inset in Fig. 2). One is the direct decay within the right 
well out of the eigenstate of the two-state Hamiltonian 
localized predominantly in the left well, but with a small 
probability amplitude in the right well. The rate of this 
decay is F (|7h. Another is the transition between the 



two eigenstates induced by the interwell relaxation that 
transfers the probability between the two wells and is 
followed by the intrawell decay out of the lower-energy 
eigenstate with the rate F2. The rate of the interwell 
transition between the eigenstates is 
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(20) 



At /i < A, the second path dominates, and for sufhciently 
large e, e » ^ I P, the bottleneck of the relaxation pro- 
cess is the interwell transition between the eigenstates 
and the tunneling rate ( p^ becomes independent of the 
intrawell rate F2: — 1- In general, the competition 
between the two types of transitions gives an expression 
for the tunneling rate, t^^ = F2(7 -I- f )/(7 -I- F2), that 
agrees with eq. (|l9| ) at e S> A. 

For large p^ p ^ the interwell relaxation is weak 
close to resonance, and for sufficiently small bias, e <C /i, 
eq. ( |l9| ) coincides with eq. (|lj) with F2 ^ A. However, 
at larger bias, e ^ p, the interwell relaxation increases 
and the tunneling rate is again given by 7 (pO[). Only 
when p becomes comparable to the level separation ujp 
in the wells, the interwell relaxation is completely negli- 
gible and the resonant tunneling peak has the Lorentzian 
shape for all relevant energies. 
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FIG. 2. The zero-temperature rate of flux tunneling 
between the wells of the double- well potential as a function of 
the energy bias e for different values of /i, which characterizes 
the relative strength of the two relaxation mechanisms (p^. 
The inset shows the diagram of the off-resonant transitions. 

Away from resonance, when e ^ Up, the interwell tun- 
neling with the rate ( pO| ) corresponds to the transitions 
between the states localized in opposite wells that are the 
closest in energy, while "intrawell" tunneling (|l^) corre- 
sponds to the transition between the states that are at 
least next-nearest neighbors in energy. Although there 
is no qualitative difference between the two types of the 
transitions away from the resonance, they lead to very 
different shapes for the tunneling peaks close to reso- 
nance. 
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IV. PHOTON-ASSISTED TUNNELING 

When the SQUID is irradiated with an external rf sig- 
nal, the macroscopic resonant flux tunneling can go more 
effectively through one of the excited states in the left 
well of the SQUID potential rather than out of the ground 
state, since the amplitude of tunneling — A/2 out of the 
excited state is much larger than the tunneling ampli- 
tude for the ground state. In this Section, we consider 
the situation when an rf signal of frequency to resonantly 
couples the ground state |0) in the left well of the poten- 
tial to an excited state |1) with energy E in this well (Fig. 
3). The energy E is on the order of Wp, and the condition 
of the resonant excitation is that the detuning v = E — lj 
is small, v <^ ujp. If the amplitude a of the rf excitation is 
also relatively small, a <C Wp, the off-resonant coupling to 
other stateg23 is not important, and the coupling between 
the states |0) and |1) can be described in the rotating- 
wave approximation. In this approximation, the terms in 
the coupling that oscillate rapidly (with frequencies on 
the order of ujp) are neglected, and the coupling Hamil- 
tonian is written as: 



states |0), and |2). Time evolution of p is described 
by the same eq. (0) but with the Hamiltonian Hq given 
now by eq. (p2|). We begin discussion of the dynamics 
of flux tunneling in the three-state system (^) with the 
case when the intcrwcll relaxation "f[p\ can be neglected, 
and the only mechanism of the energy relaxation in the 
system is the intrawell relaxation T[p\. This relaxation 
is characterized by the two rates Ti^2 of the transitions 
from the states |1, 2) into the lower energy states of the 
left and right potential well, respectively (Fig. 3). An 
obvious generalization of eq. (^) for T[p] to the three- 
state basis gives the off-diagonal part of eq. ^ for the 
time evolution of p of the following form: 

Poi = («i^ - ri/2)poi + iaipoo - pii)/2 - iApo2/2 , 

Pi2 = -{is + (Fi + F2)/2)pi2 + iA(p22 - Pii)/2 - mpo2/2 , 

P02 =- {iiv -£)- F2/2)po2 - WP12/2 - 2Apoi/2 . (23) 

Equations (|2^) allow us to express the off-diagonal ele- 
ments of p in terms of the diagonal ones in the stationary 
regime. Inserting the stationary values of the off-diagonal 
elements into the equations for the diagonal elements: 



|l)(0|e+*'^'). 



(21) 



If the excited state |1) in the left well is coupled reso- 
nantly with amplitude —A/2 to a state |2) in the right 
well that is shifted in energy by e with respect to |1), the 
total Hamiltonian for the flux dynamics in the basis of 
the three states |0), |1), and |2) is: 
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a/2 


V 









(22) 
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FIG. 3. Schematic diagram of the "photon-assisted" 
macroscopic resonant tunneling of flux stimulated by an rf 
perturbation of strength a/2. 

As in the previous Section, the average flux tunneling 



rate r can be calculated according to eq. (^ from the 
stationary density matrix p of the system in the basis of 



Poo = -almpoi + Fipii -|- F2P22 , 
Pll = almpoi + AImpi2 - Fipn 
P22 = -AImpi2 - F2P22 , 



(24) 



we calculate the stationary probability P22^ and find the 
fiux tunneling rate (|l5|). 

Equations ( p^ were written under the assumption that 
the relaxation in the left well brings the system out of the 
state |1) directly into the ground state |0). Although this 
is strictly true only in the case when 1 1) is the first excited 
state in the well, eqs. (2^) can be also used to calculate 
the average fiux tunneling rate in other situations. In- 
deed, when the rf signal drives the system into the state 
|1) that is not the first excited state, the intermediate 
states in the left well that exist between the states |0) 
and |1) are populated by the process of relaxation out of 
|1). If the flux tunneling out of these states is neglected 
(similarly to tunneling out of the state |0)), their effect 
on the average tunneling rate can be accounted for by in- 
clusion of the occupation probabilities of these states in 
the normalization condition. Since these probabilities in 
the stationary regime are proportional to the stationary 
occupation probability pn, this can be done through an 
additional factor A in the normalization condition for the 
state |0), |1), and |2): 



Poo + Apii + P22 = 1 ■ 



(25) 



To give an example, we can calculate the factor A assum- 
ing that the left well is parabolic in the relevant energy 
range. Then, the standard result for the-|linear relax- 
ation of the harmonic oscillator (see, e.g.,E^) is that the 
oscillator makes the transitions only between the nearest- 
neighbor states and that the transition rate from the 
state |m) into |m — 1) is proportional to to. This means 
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that in the stationary state pm-i,m-i = fnpm,m/ {m— 1), 
and if the rf radiation drives the system into the nth ex- 
cited state of the left well, then A = "-Sm=i(V"^)- To 
avoid extra parameters, however, we assume from now 
on that the state |1) is the first excited state in the left 
well, so that A = 1. 

Sufficiently compact analytical expressions for the tun- 
neling rate can be obtained from eqs. ( |2^ ) and (|2j) only 
in certain limits. For example, for smalLtf amplitude and 
weak relaxation, a ^ ri.2 ^ A, we gettj: 

{2i^-e- f7)2(2!y - e + fl)^ + 4(ri(i^ - e) + ra//)^ ' 

(26) 

Equation (^) describes two peaks (discussed in more de- 
tails below) in the dependence of the tunneling rate on 
the detuning v. The peaks are broadened by the relax- 
ation, and their positions correspond to the two eigen- 
states of the Hamiltonian (|^): v ^ {e ± 51)/2. Another 
expression for the tunneling rate can be obtained for large 
relaxation rates Fi 2 ^ a. A: 



(4j.2+r2)(4(i.-£)2+F2) ■ ^ ' 

Depending on the relation between the energy bias e and 
the relaxation rates, the tunneling rate (|2^) as a function 
of detuning v contains either one (for small e) or two 
separate peaks (for large e). The peak positions in this 
case coincide with the position of the energy levels 1 1 , 2) 
localized in the two wells. 

For arbitrary parameters it is convenient to use the sta- 
tionary solution of eqs. ( ^ ) and (|2j) to plot the tunneling 
rate numerically. To simplify the discussion, we as- 
sume that the relaxation rates in the two wells are the 
same, Fi = F2 = F. Figure 4 shows the dependence of 
on detuning v obtained in this way for several rf am- 
plitudes a in the regime of small relaxation rate F. The 
main qualitative feature of Fig. 4 (that can also be seen 
in eq. (|26|)) is that for small a the resonant peak in the 
tunneling rate r is split into two peaks due to coherent 
oscillations of flux between the two wells. The appear- 
ance of such splitting can be easily understood, since a 
weak rf signal excites the system not into the state |1) lo- 
calized in the left well, but into the two hybridized states 
formed out of the states |1) and |2) in the two wells. The 
two different energies of the two hybridized states lead 
to the two peaks in the tunneling rate. This means that 
the splitting of the resonant-tunneling peak is the direct 
manifestation of the quantum coherent oscillations of flux 
between the two wells of the SQUID potential. 

Figure 4 illustrates also how the splitting of the res- 
onant tunneling peak is suppressed with the increase of 
the rf amplitude. Suppression of the peak splitting can 
be understood in terms of the time 1/A required to es- 
tablish the stationary states hybridized between the two 
wells. Large rf amplitude a on the order of the tunnel 



amplitude A causes rapid Rabi oscillation between the 
states |0) and |1) and does not allow for sufficient time 
to establish the two hybridized states. The system is 
therefore effectively excited into the state |1) localized 
in the left well, and single resonant tunneling peak is 
formed around the energy of this state. Broadening and 
suppression on this single peak seen in Fig. 4 as a in- 
creases beyond A, is a version of the generic "quantum 
Zeno" effect, when tunneling out of a metastable state is 
suppressed by rapid perturbation of this state. 



^ 0.2 



a/A=LO 




r=o.iA 


f\/\0.6 


e=0.0 


lo.sl 


2.0 i 




3-0 Jfn 







-4.0 -2.0 0.0 2.0 4.0 



v/A 

FIG. 4. The rate of the photon-assisted resonant tun- 
neling of flux $ between the two weUs of the SQUID potential 
as a function of detuning v for several values of the rf ampli- 
tude a. For discussion see the text. 

Figure 5 shows the evolution of the coherently-split 
tunneling peaks at small rf amplitude a with the bias en- 
ergy £ and with the relaxation rate F. We see that with 
increasing energy bias (Fig. 5a) , the peaks follow the po- 
sition of the energy levels and the splitting between them 
increases. Simultaneously, the peak height decreases re- 
flecting the overall suppression of the tunneling rate as 
one moves away from the resonance. Figure 5b shows how 
the double-peak structure in the small-bias regime repre- 
senting the coherent mixing of the flux states in the two 
wells is suppressed by increasing relaxation rate F. The 
structure is visible up to the relaxation rates F ~ 0.5A. 

The peaks in Fig. 5 are shown only for e > 0. The 
peak structure for negative e can be understood from the 
"symmetry" relation t~^(— t^, — e) = r^^iy^e) that can 
be deduced from eqs. ( p3| ) and (p^). Equations (|2^ ) show 
that in the stationary regime, changing the sign of e is 
equivalent to changing the sign of a, A and replacing the 
off-diagonal elements of p with their complex conjugate 
values. This transformation obviously does not change 
the transition rates in eq. (|2^), and therefore does not 
change the flux tunneling rate r^^. 

Another interesting manifestation of the coherent flux 
tunneling between the wells can be seen in the depen- 
dence of the tunneling rate on the bias energy e at fixed 
detuning v (Fig. 6). For weak relaxation, the hybridized 
states are well-developed, and when the rf excitation en- 
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ergy lies between these two states, the tunneling rate is 
strongly suppressed. At = 0, when the system is ex- 
cited precisely into the state |1) localized in the left well, 
this condition is satisfied and tunneling rate is strongly 
suppressed for any energy bias e. As can be seen from eq. 
( Pq ) and Fig. 6, in this case the tunneling rate as a func- 
tion of e is described by a Lorcntzian centered around 
e = 0. In contrast to resonant tunneling peaks in the in- 
dependence of the tunneling rate, which have small width 
proportional to the relaxation rate F, the width of this 
Forcntzian is large, A^/F, and is inversely proportional 
to F. When the detuning v deviates from zero, there 
is an energy bias e at which the excitation energy co- 
incides with the energy of one of the hybridized states. 
The tunneling rate has a peak under such resonance con- 
ditions. For not-too-small i^'s, this resonant peak again 
has a small width proportional to F. In Fig. 6, one can 
see how the transition between the broad and narrow 
tunneling peaks takes place for v > Q. As before, the 
results for negative detuning can be deduced from the 
relation T~-'^(—;/, —e) =T~^{y^£). 
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FIG. 5. Evolution of the double-peak structure in the rate 
of the photon-assisted resonant flux tunneling as a func- 
tion of detuning u with increasing (a) bias energy e, and (b) 
relaxation rate F. 

Finally, we discuss the effect of weak interwell relax- 
ation on the photon-assisted tunneling. We start by gen- 



eralizing eq. (^) for this relaxation to the three-state sit- 
uation relevant for the photon-assisted tunneling. Under 
the natural assumption that the average flux in the states 
|0) and |1) in the left well of the SQUID potential is the 
same, the part of the dissipative coupling (||) that corre- 
sponds to the interwell relaxation is: 



V = -IfS<^ 




(28) 



While eq. (|2|) is written in the flux basis |0) , |1) , |2), 
weak relaxation is conveniently described in the basis of 
the eigenstates \n) of the Hamiltonian (p2|). In this basis, 
the contribution of the interwell relaxation ( p8| ) to the 
evolution of the density matrix p is given by the standard 
expression similar to eq. (^: 

Pnn — ^ ^^ {pininPmm ^nmPnn) j (^'^) 
m 

Pnni = - Wmn + ^ ^('^^ik + lrnk)]Pnm , U ^ m . 
k 

Transition and dephasing rates in these equations are: 



In 



1 



-(e„-e™)/T ' 



In 



where Unm are the matrix elements of the operator U 
( p8| ) in the eigenstates basis, and £„ is the energy of the 
eigenstate \n). 
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FIG. 6. Bias-energy dependence of the photon-assisted flux 
tunnehng rate at fixed detuning v for small rf amplitude and 
relaxation rate F. For vanishing detuning, the tunneling rate 
exhibits a very broad maximum with a width inversely pro- 
portional to F. 

Interwell relaxation can be included into the evolution 
equations for the density matrix on the basis of eq. ( [29| ) 
numerically. We diagonalize the Hamiltonian (|2^), calcu- 
late the interwell relaxation terms (|2^) in the eigenstates 
basis, and transfer them into the flux basis, where the 
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intrawell relaxation has the simple form (p3|), (24). Cal- 
culating finally the stationary value of the density matrix 
p we find the flux tunneling rate (15). 

Figures 7 and 8 show results of such a calculation ob- 
tained at vanishing temperature T. In Figure 7, the tun- 
neling rate is plotted as a function of the detuning v for 
e = and several values of the relative strength of the 
intrawell relaxation F. The eigenstates of the two-state 
Hamiltonian at e = are symmetric between the two 
wells, and in absence of the interwell relaxation produce 
two symmetric resonant tunneling peaks (see Figs. 4 and 
5a). As can be seen from Fig. 7, the interwell relaxation 
makes the tunneling peaks asymmetric. The positive-;^ 
side of the double-peak structure corresponds to excita- 
tion of the system into the lower-energy eigenstate and 
is unaffected by the interwell relaxation at zero temper- 
ature, since there is no energy in this regime to create 
additional tunneling path. In contrast, the negative-;^ 
side of the double-peak structure corresponds to excita- 
tion of the system into the eigenstate with larger energy, 
and the interwell relaxation increases the rate of tunnel- 
ing out of this state. Because of this, the negative-i' peak 
in Fig. 7 is larger than the tunneling peak at positive 
and the tunneling rate at < decreases much more 
slowly away from the peak than at positive v. 
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FIG. 7. The rate t^^ of the photon-assisted flux tunnehng 
as a function of the detuning in the case of symmetric cou- 
pling between the tunneling flux states, e = 0, in the presence 
of both interwell and intrawell relaxation. Different curves 
correspond to different magnitudes of the intrawell relaxation 
rate F relative to the interwell relaxation rate. 

Interwell relaxation introduces asymmetry also in the 
dependence of the flux tunneling rate on the bias energy 
£. Examples of such dependence are shown in Fig. 8 for 
two values of the detuning, u = ±2A. In both cases, 
the tunneling rate has a resonant peak at e ~ i/ similar 
to the peaks shown in Fig. 6 for vanishing interwell re- 
laxation, when the peaks for the two values of detuning 
are symmetric. Comparison of Figs. 8a and 8b shows 
that the interwell relaxation makes the peaks asymmet- 
ric. In particular, the peak at < is smaller than the 



peak at 1/ > 0. Although this asymmetry appears to be 
opposite to that in Fig. 7, where the negative-i^ peak is 
larger, it has the same origin as in Fig. 7. The peaks at 
V < Q and v > Q correspond to resonant excitation of 
the system into, respectively, the upper and lower energy 
eigenstates. When the resonance occurs for |e| > A, as 
in Fig. 8, the eigenstates are already to a large extent 
localized in one or the other well. At e ~ v > 0, the 
lower eigenstate is centered in the right well and the in- 
terwell relaxation increases the tunneling rate, while at 
e ~ 1/ < the lower eigenstate is centered in the left well 
and the interwell relaxation brings the system back to 
this well suppressing the tunneling rate. As a result, the 
resonant tunneling peak in Fig. 8a {v > 0) is larger than 
in Fig. 8b {v < 0). The height of the negative-i/ peak is 
more sensitive to the relative strength of the two relax- 
ation mechanisms and decreases with decreasing rate F 
of the intrawell relaxation. 
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FIG. 8. The rate of the photon-assisted flux tunneling 
as a function of the bias energy e for two "symmetric" values 
of the detuning f and the same rates of the intrawell relax- 
ation r as in Fig. 7. In (a), decreases with increasing F 
for e above the peak, while in (b), the peak height increases 
with F. 

The tails of the photon-assisted resonant peaks can 
also be described analytically. When the bias energy and 
detuning are not close to any resonance, \i>\, |e|, |!^ — £| S> 
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a, A both the interwell tunnehng and rf excitation can be 
treated as perturbations. The dynamics of flux tunnehng 
in this regime can be described as a coexistence of the 
two tunnehng paths similar to the off-resonant tunnehng 
discussed in Sec. 2 (see inset to Fig. 2). li i' > e, the 
effective energy v — e the state |2) in the right weU 
(i.e., the energy of this state brought down by a quantum 
of rf radiation, as in the Hamihonian (^)) is above the 
energy of the initial state |0) in the left well, and the only 
energy- allowed tunneling path is direct relaxation in the 
right well out of the perturbed state |0). Similarly to eq. 
(O), perturbation theory in a, A gives for the rate F of 
this tunneling: 



F^g^A^ 
lQv'^{v — e 



(30) 



If < e, the effective energy of the state |2) is lower 
than the energy of the state |0), and in addition to the 
tunneling ( pO[ ) there is a competing tunneling process. It 
consists of a transition between the perturbed states |0) 
and 1 2) with the rate 7, 



7 ^ 



(31) 



that is driven by interwell relaxation, followed by direct 
relaxation in the right well with the rate F2. As in Sec. 
2, the coexistence of the two tunneling paths gives the 
following total tunneling rate: 



r =12: 



"7 + r2 



(32) 



One can check that eqs. ( |30| ) and (^) agree, respectively, 
with the negative-e tail of the resonant tunneling peak in 
Fig. 8b and the positive-e tail of the peak in Fig. 8a. 

In summary, we have studied the effects of two types 
of relaxation mechanisms on the macroscopic resonant 
tunneling of flux in SQUIDs under stationary-bias condi- 
tions and with external rf irradiation. Coherent splitting 
of the resonant-tunneling peaks by rf radiation provides 
a convenient way of studying quantum coherence of flux 
states. 

This work was supported by ARO grant 
DAAD199910341. 

Appendix 

In the Appendix we show explicitly how the Hamilto- 
nian of the two-well system (n|) can be reduced at reso- 
nance to the two-state form (0), and derive an expression 
for the tunneling amplitude A. Assuming that the trans- 
parency D of the barrier separating the two wells is small, 
D <C 1, and that the resonance occurs between the states 
with large n, we can use the WKB approximation for the 
wavefunctions '0jn($) of the Hamiltonian (^. In this ap- 
proximation, the wavefunction between the right and the 
left turning points r and / is: 



■^{^) = —= cos(w(<I>) - (5) , 



(33) 



with the WKB phase w($) = (1/^) /* d$'p(^') - 7r/4 
and momentum p = \1C{E - f7($))]i/2. in eq. (|3|), b 
is a constant phase shift, and ^ is a normalization con- 
stant. The phase w is defined in such a way that (5 = for 
an isolated well, when the wavefunction decays exponen- 
tially in the classically inaccessible region $ > r, $ < L 
When the two wells are coupled, the energy E of the 
state common to them deviates from the eigenenergies 
En of the isolated wells that are determined by the Bohr- 
Sommerfeld condition (l/?i) J[ pd^ = 27r(n -I- 1/2). For 
weak tunneling, this deviation is small in comparison to 
the state energy, and creates a small but non-vanishing 
phase shift 6: 



6=1 f {p{E) - p{e,,)) = n— 



(34) 



where ujp is the frequency of the classical oscillations in 
the well and which in the WKB approximation deter- 
mines the spacing of the energy levels. 

At non-vanishing 5, the wavefunction (|3^ ) has a part 
that grows exponentially in the classically inaccessible 
region, as one can see rewriting it as 



Ae 



iS 



VP 



cosw{^) 



iA sin S 



,iu)(<I>) 



(35) 



According lo general rules of the WKB approximation 
- see, e.g.,LEl, the two terms in this expression produce 
exponentially decaying and growing components of the 
wavefunctions with the amplitudes Ae^^ /2 and —AsinS, 
respectively. To find the energy E we need to match the 
amplitudes of the wavefunctions of the right and left wells 
in the barrier region. Under the conditions of resonance, 
equating the amplitudes and keeping only the terms of 
the first order in D we get: 

A,^^ = -A,^, A,^^ = ^A,^, (36) 
LOi In UJ2 27r 

where e^ , j = 1, 2 are the energies of the resonant states, 

andD = eM-i^/^)Irl\p\d^- 

The probability to be in the right/left well is directly 
related to the wavefunction amplitudes Aj: 



B\A, 



(37) 



where B is the j-independent part of the normalization 
constant. Introducing the amplitudes a of this probabil- 
ity, aj = Aj/ujy^, we bring eq. ( p6| ) into the form that 
coincides with the Schrodinger equation of the two-state 
system: 

Eai = e\a\ — — a2 , Ea-i = £2(^2 — —ai , (38) 
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with the tunnehng amplitude 



A 



1/2 



-D. 



Since the resonant states |1, 2) are orthogonal to all other 
states of the full Hamiltonian this proves that their 
dynamics can be described by the Hamiltonian (|^). 

Away from resonance, when e = ei — £2 ^ ujp, one 
of the amplitudes a (for the states localized in the left 
well, 02) is small, a2 ~ D. Keeping, as before, only 
the terms of the first order in D, we see from eq. ( |3q ) 
that in this case Ei = ei, i.e. Si = 0. Matching the 
growing wavefunction in the right well with the decaying 
wavefunction in the left well we find that away from the 
resonance 



ttA 



Q!2 



2uj2 sin(7re/a;2) 



(39) 



Equation (^) extrapolates smoothly between the succes- 
sive resonances. 

Making use of the harmonic approximation for the po- 
tential [/(<&), one can show that the results of this Ap- 
pendix (eqs. (38) and (|39|)) can be extended to the low- 
lying states with small n which cannot be described with 
the WKB approximation. Such an extension leads only 
to minor modifications in the definition of the barrier 
transparency D in these equations. 
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